#ifndef BOSS_Afterburner_ReconstructedParticle_H
#define BOSS_Afterburner_ReconstructedParticle_H

#include "Fit/FitPars.h"
#include "Particle.h"
#include "RooDataHist.h"
#include "RooGaussian.h"
#include "RooRealVar.h"
#include "TString.h"
#include <list>
#include <utility>
#include <vector>
#include <yaml-cpp/yaml.h>

/// @addtogroup Afterburner_objects
/// @{

/// Extended data container of a reconstructed particle.
/// In essence an container of `ROOT`'s `TParticlePDG`, but with extended functionality for particle reconstruction. It contains additional parameters that are relevant when reconstructing the particle, like fit range in the invariant mass plot. The idea is that you can use this object to generalise functions like fitting, where the mass range is particle-specific. It is then convenient to only have to feed the function this object.
/// @todo Somehow, a distinction has to be made for different decay channels (not just daughters) in the parameters of `ReconstructedParticle`.
/// @author   Remco de Boer 雷穆克 (r.e.deboer@students.uu.nl or remco.de.boer@ihep.ac.cn)
/// @date     October 25th, 2018
class ReconstructedParticle : public Particle
{
public:
  /// @name Constructors
  ///@{
  ReconstructedParticle() {}
  ReconstructedParticle(Int_t pdgCode);
  ReconstructedParticle(const Char_t* pdgCode);
  ReconstructedParticle(Int_t pdgCode, const Char_t* daughters);
  ReconstructedParticle(const Char_t* pdgCode, const Char_t* daughters);
  ReconstructedParticle(const YAML::Node& input);
  ///@}

  /// @name Setters
  ///@{
  void SetDaughterLabel(const Char_t* daughters);
  void SetDaughterLabel(const YAML::Node& input);
  void SetFitPars(const TString& input_sig, const TString& input_bkg)
  {
    fFitParameters.Set(input_sig, input_bkg);
  }
  void SetGaussianPars(TString input);
  void ImportFitPars_DG();
  void SetFitRange(Double_t from, Double_t to) { fFitRange = {from, to}; }
  void SetPlotRange(Double_t from, Double_t to) { fPlotRange = {from, to}; }
  void SetFitRange(std::pair<Double_t, Double_t> range) { fFitRange = range; }
  void SetPlotRange(std::pair<Double_t, Double_t> range) { fPlotRange = range; }
  ///@}

  /// @name RooFit
  ///@{
  std::shared_ptr<RooRealVar>              RooRealVarInvMass();
  std::shared_ptr<RooRealVar>              RooRealVarMean();
  std::shared_ptr<RooRealVar>              RooRealVarMeanZero();
  std::shared_ptr<RooRealVar>              RooRealVarSigma(UChar_t num);
  std::shared_ptr<RooRealVar>              RooRealVarSigmaFixed(UChar_t num);
  std::vector<std::shared_ptr<RooRealVar>> FitConvolutionBWGaussian(TH1F* hist, TString logScale);
  std::vector<std::shared_ptr<RooRealVar>> FitBreitWigner(TH1F* hist, TString logScale);
  std::vector<std::shared_ptr<RooRealVar>> FitPureGaussians(TH1F* hist, TString logScale,
                                                            TString    outputName = "",
                                                            const bool fixpars    = false);
  RooDataHist CreateRooFitInvMassDistr(TH1F* hist, const RooRealVar& var);
  ///@}

  /// @name Getters
  ///@{
  const Char_t*  DaughterLabel() const { return fDaughterLabels.Data(); }
  const Double_t BWConvolutedWidth() const { return fBWConvolutedWidth; }
  const Double_t BWPureWidth() const { return fBWPureWidth; }
  const Double_t FitFrom() const { return fFitRange.first; }
  const Double_t FitUntil() const { return fFitRange.second; }
  const Double_t GaussianSmallWidth() const { return fGaussianWidths.front(); }
  const Double_t GaussianWideWidth() const { return fGaussianWidths.back(); }
  const Double_t GaussianWidth(size_t i) const { return fGaussianWidths.at(i); }
  const Double_t LowerMass() const;
  const Double_t MassOffsetPercentage() const { return fMassOffset; }
  const Double_t PlotFrom() const { return fPlotRange.first; }
  const Double_t PlotUntil() const { return fPlotRange.second; }
  const Double_t UpperMass() const;
  const Char_t   NPol() const { return fNPol; }
  const FitPars& FitParameters() const { return fFitParameters; }

  const std::vector<Double_t> GaussianWidths() const { return fGaussianWidths; }
  ///@}

protected:
  /// @name Data members
  ///@{
  FitPars fFitParameters;
  ///< Vector of fit parameters. @remark Temporary implementations, see `FitPars`.
  TString  fDaughterLabels; ///< LaTeX formatted string for decay particles.
  Double_t fBWConvolutedWidth;
  ///< Estimate for the width of the Breit-Wigner function when convoluted with a Double_t Gaussian.
  Double_t fBWPureWidth;
  ///< Estimate for the width of the Breit-Wigner function when fitting BW only.
  Double_t fMassOffset; ///< Percentage (divided by 100) that the mean (namely, the mass) may vary.
  Char_t   fNPol; ///< Degree of the background polynomial (whether Chebychev or normal polynomial).
  std::vector<Double_t> fGaussianWidths;
  ///< Vector of for Gaussian widths that you can use that as estimates of the widths for e.g. the
  ///< Double_t gaussian that you plan to fit. These sigmas are supposed to characterise the
  ///< resolution of the detector. For consistency in naming, the first one should be smaller than
  ///< the second and so forth.
  std::pair<Double_t, Double_t> fFitRange;
  ///< Invariant mass range over which you fit a function (Double_t Gaussian, Crystal ball,
  ///< Breit-Wigner, etc.).
  std::pair<Double_t, Double_t> fPlotRange;
  ///< Invariant mass range that you plot.
  std::list<std::shared_ptr<RooRealVar>> fRooRealVars;
  ///< This list is to ensure that `RooRealVar` remain existing untill this `ReconstructedParcticle`
  ///< is destroyed.
  std::list<std::shared_ptr<RooGaussian>> fRooGaussians;
  ///< This list is to ensure that `RooGaussian` remain existing untill this
  ///< `ReconstructedParcticle` is destroyed.
  ///@}

  /// @name Helper methods
  ///@{
  void DetermineReconstructionParameters();
  void Initialize();
  ///@}
};

/// @}
#endif